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ABSTRACT 

We present the detailed global structure of black hole accretion flows and outflows through newly 
performed two-dimensional radiation-magnetohydrodynamic simulations. By starting from a torus 
threaded with weak toroidal magnetic flelds and by controlling the central density of the initial torus, 
po, we can reproduce three distinct modes of accretion flow. In model A with the highest central 
density, an optically and geometrically thick supercritical accretion disk is created. The radiation 
force greatly exceeds the gravity above the disk surface, thereby driving a strong outflow (or jet). 
Because of the mild beaming, the apparent (isotropic) photon luminosity is ~ 22Le (where Le is 
the Eddington luminosity) in the face-on view. Even higher apparent luminosity is feasible if we 
increase the flow density. In model B with a moderate density, radiative cooling of the accretion 
flow is so efficient that a standard- type, cold, and geometrically thin disk is formed at radii greater 
than ~ 7Rs (where Rs is the Schwarzschild radius), while the flow is radiatively inefficient otherwise. 
The magnetic-pressure-driven disk wind appears in this model. In model C the density is too low 
for the flow to be radiatively efficient. The flow thus becomes radiatively inefficient accretion flow, 
which is geometrically thick and optically thin. The magnetic-pressure force, in cooperation with the 
gas-pressure force, drives outflows from the disk surface, and the flow releases its energy via jets rather 
than via radiation. Observational implications are briefly discussed. 

Subject headings: accretion, accretion disks - black hole physics - ISM: jets and outflows - magneto- 
hydrodynamics (MHD), - radiative transfer 
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1. INTRODUCTION 

Black hole accretion disks provide the most power- 
ful energy-production mechanism in the universe. How- 
ever, theoretical development in this area is rather be- 
hind, compared with that of the stars. The central en- 
ine of the dis ks was identifled as having magnetic origin 
Balbus fc Ha wlev 1991) in the late 1990s after extensive 
investigations had been carried out by many authors. 
Yet, no complete picture of magnetized accretion flow 
and outflow has yet been obtained at this moment. We 
may thus say that the theory of the accretion disks is 
now in a similar situation to that of the stars in the 
1940s-1950s. 

Before the identiflcation of the central engine, one- 
dimensional accretion disk models were constructed 
based on the phenomenological a- viscosity prescription, 
whereby the viscous torque is proport ional to the pres- 
sure, after the pioneering work by (Shakura & Sunyaev| 
IT973). The standard disk model was flrst established as 
a model for accretion disks at moderately high luminosi- 
ties, followed by various disk models, including the slim 
disk model and the radiatively inefficient accretion flow 
(RIAF) model, which are proposed as models for the ac- 
cretion flow with hig:her Le) and lower luminosities , 
respec tively (Shakura & Sunvaev 1973; Ichimaru' '1977"; 



.ees e t al. 1982; Abramowicz et al. 1988; Narayan & Yi 
1994 see iKato et al.l l2QQ8l for an extensive review). 



the observations (e.g., through spectral fltting), how- 
ever, the results derived those simplifled models need 
to be checked, since their results may depend on the 
a- viscosity assumption and radially one-dimensional ap- 
proximation. In fact, Hirose et al.l (2009) have claimed 
that the radiation-pressure-dominated part of the disk is 
thermally stable, though it was shown to be unstable 
unde r the a-viscosity pre scription (;Shibazaki fc Hoshil 
[T975I : iShakura fc Sunvaevl [T976). 

As a distinct line of disk research, multi- 
dimensional hydrodynamical flow simulations were 
attempted based on the a-viscosity prescription 
jlgumenshchev fc Abramowiczl 119991 : iStone et all 119991 : 
McKi nnev fc Gammid [20021 ). but it was not a main 
stream of research. After the identiflcation of the disk 
viscosity, the situations drastically changed. Global 
magnetohydrodynamic (MHD) simulations have be en ex- 
tensive ly performed by several groups (Matsumot^Hl ' 



These models have good prediction powers and thus 
make it possible to directly compare the theory with 



Machid a et al. 2000; Hawlev & Krolik 2001; Koide e t^ 
20011 : iDe Villi ers et al. 2003, Hawlev & Krolil^ .2006 
see, however, a review by Spruit 2010). Despite these 
studies, there is no wide consensus regarding the global 
and local behavior of magnetic flelds; e.g., it is an 
open question how numerical results are sensitive to 
the initial magnetic fleld strengths and conflgurations, 
boundary conditions, and numerical resolutions. This is 
due partly to highly nonlinear spatio-temporal evolution 
of magnetic flelds in differentially rotating media. 
Further, most of the previous global MHD simulations 
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are non-radiative ones and, hence, they cannot model 
accretion disks in high- luminosity states, in which 
significant radiative cooling (and occasionally strong 
matter-radiation coupling) is expected. The radiation 
transfer should be solved to explain the energy release 
processes within the disk. The dynamical effects of 
radiation are especially important for the radiation- 
pressure-dominated disk, since they are expected to 
produce strong radiatively driven outfiows, thereby the 
disk structur e being significantly modified. 

In Ohsuga et_alj ([2009) we have presented a new 
type of accretion fiow simulation based on the global 
radiation-magnetohydrodynamic (RMHD) simulations. 
We solved the problem from the first principle (i.e., 
without employing the phenomenological a- viscosity pre- 
scription), including the case of radiatively very efficient 
accretion fiow. That is, we considered the following ba- 
sic processes in the accretion flow and jets: the trans- 
port of angular momentum induced via the magnetic 
torque, leading to the accreting motion, the conversion 
of the mechanical energy to the thermal energy via the 
MHD processes, the dissipation of thermal energy, the 
radiative transfer, and radiation-pressure and Lorentz 
forces, which play important roles in launching outflows 
and supporting the disks in the vertical direction. The 
overview of the RMHD simulation results of the accretion 
flow and outflow around the black holes has already been 
published elsewhere (e.g., Ohsuga et al. 2009). Here, we 
present the detailed analyses of the simulated flow struc- 
ture based on our newly performed, improved simula- 
tions. 

Similar RMHD simulations were attempted previously 
but only under the shearing-box approxima tion, in which 
only a local patch of the disk is treated (e.g. jTurner et al.] 
120031 : ^ose et al. 2006). Global (radial) coupling of 
magnetic flelds are ignored in those simulations and thus 
collimated outflows cannot be produced there. Further, 
advective motion of gas and photons were not considered. 
For establishing a realistic picture of accretion disks and 
outflow, therefore, global, multi-dimensional RMHD sim- 
ulations, of a kind reported in the present paper, are in- 
dispensable. We will show that the realistic flow proper- 
ties simulated here share some similarities with the flows 
described by the previous one-dimensional models but 
that they also exhibit new features which were not an- 
ticipated previously. The plan of this paper is as follows. 
Basic equations and physical assumptions are explained 
in Section 2, and numerical procedures are described in 
Section 3. We will then show the results of simulations 
in Section 4 and give discussion in Section 5. Finally, 
Section 6 is devoted to summary. 

2. BASIC EQUATIONS AND ASSUMPTIONS 

We solve a full set of RMHD equations under flux- 
limited diffusion (FLD) approximation in the cylindrical 
coordinates, (r, (/p,z). In the present study, we assume 
that the flow is non-self-gravitating, reflection symmetric 
relative to the equatorial plane (with z = 0), and axisym- 
metric with respect to the rotation axis (i.e., d/d^p = 0). 
We describe the gravitational fleld of the black hole in 
terms of pseudo-Newtonian hydrodynamics, in which the 
gravitational potential is given hy t/j = —GM/{R — Rs) 
(jPaczvnskv &: Wiital [T980h . where R[= (r^ + z^)^/^] is 
the distance from the origin and Rs{= 2GM/c^) is the 



Schwarzschild radius (with M and c being the black hole 
mass and the light velocity, respectively). 
The basic equations are the continuity equation. 



dp 

dt 



V • {pv) = 0, 



(1) 



the equations of motion, 

d{pv) ^ ( BB 
^ ' V ■ I pvv 



dt 



-V p 



Stt 



-^Fo-pVi^, (2) 

c 



the energy equation of the gas, 

— + V ■ {ev) = -pV ■ V + -^VJ - 47r/tB + ckEq, (3) 
the energy equation of the radiation, 

dt 



V • (Eov) = - V • Fo - Vv : Po + ^nnB - ckEq, 

(4) 



and the induction equation. 



— — = V X { V X B 77 J 

dt V c ' 



(5) 



Here, p is the gas mass density, v = {vr^Vi^p^Vz) is the 
velocity, e is the internal energy density of the gas, p is 
the gas pressure, B = {Br, B^^B^) is the magnetic fleld, 
J(= cV X B/Ati) is the electric current, r] is the resis- 
tivity, B is the blackbody intensity, is the radiation 
energy density, Fq is the radiation flux, Pq is the radia- 
tion pressure tensor, n is the absorption opacity, and x 
is the total opacity. 

The gas pressure is related to the internal energy den- 
sity of the gas by 



P= (7- 



(6) 



where 7 is the speciflc heat ratio. The temperature of 
the gas, Tgas, can then be calculated from 



P 



pkT, 



gas 



/i777p 



(7) 



where k is the Boltzmann constant, fi is the mean molec- 
ular weight, and 777p is the proton mass. 

For the absorption opacity, we consider the Rosseland 
mean free— free absorption, , and bound— free absorp- 
tion for solar metallicity, A^bf, 

= A^ff + /^bf , (8) 
where and are given by 

Kff = 1.7 X ( ] cm-1 



(|Rvbicki fc LightmanllT979h . and 



Kbf = 4.8 X 10-24r-7/2 ( 



(9) 



(10) 



(jHavashi et al.|[T962[ ). The total opacity is given by 



(11) 
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with (Jt being the Thomson scattering cross-section. 

We en iploy the FLD app ro ximat ion devel- 
oped by | Le vermore fc PomraningI (flOSTI ) ( see also 
iTurner fc Sto ne 2001). This assumption is valid in both 
the optically thick diffusion limit and the optically thin 
free-streaming limit in one-dimensional space. In this 
framework, the radiation flux is expressed in terms of 
the gradient of the radiation energy density via 



Fo = -—\/Eo. 

X 



(12) 



Here, the dimensionless function (A), which is called the 
flux limiter, is given by 



A = 



2 + 7^ 



6 + 37^ + 7^2 ' 



(13) 



using the dimensionless quantity, 1Z = | V£^o| / (x^o)- 
The radiation pressure tensor is written as 

Po = f£^o. (14) 
Here, f is the Eddington tensor and its components are 

f =i(l-/)I+i(3/-l)nn, (15) 

where / is the Eddington factor, 

/ = A + A^7^^ (16) 

and n is the unit vector in the direction of the radiation 
energy density gradient. 



iV^ol 



(17) 



We find A 1/3 and / 1/3 due to 7^ ^ in the 
optically thick limit. On the other hand, in the optically 
thin limit of 7^ ^ oo, we have |Fo| = cEq. These give 
correct relations in the optically thick diffusion limit and 
the optically thin streaming limit, respectively. 

Throughout the present study, we assume M = lOM©, 
7 = 5/3, and ja = 0.5. The resistivity is assumed to be 
constant, r] = 10~^cRs, although we have employed the 
anomalous resist ivitv (Yokovama fc Shibat a 1994) in the 
previous paper (jOhsuga et al.ll2009f ). Note that our re- 
sults do not change so much by switching the resistivity. 
The uniform resistivity might have an advantage over 
the anomalous resistivity in longer simulations, since it 
works to dissipate the magnetic fields and prevent the 
Alfven speed from being too large. Two-dimensional 
non-radiative MHD simulatio ns using a uniform resis- 
tivity have been performed by iKato et aTl ()2004a| ) . They 
succeeded in reproducing the semirelativistic jets from 
the accretion disks. 

Here we stress that in our RMHD simulations the mag- 
netic torque is responsible for the angular-momentum 
transfer and the Joule heating for the heating of 
the matter. In the conventional disk models and 
in the radiation hydrodynamic simulations, in con- 
trast, the phenomenological viscosity (a- viscosity) in- 



energy 


f dissipation ( 


Effsum et all 119881: 


Okuda fc Fuiital 


|200Ct 


Ohsusa et al.l 


2005bl: IOhsuffall2006 


). The radiative 



cooling and radiation-pressure force are both considered 
in the present study, while they cannot be taken into 
account in the non-radiative MHD simulations. 



3. NUMERICAL METHODS AND MODELS 
3.1. Outline 

Before presenting detailed calculation methods we out- 
line the procedure of our simulations in this subsection. 
We solve the time evolutions of a torus by solving the 
basic radiation-MHD equations under the assumption of 
the FLD (see Section 2). The initial torus is in hydro- 
static balance and is surrounded by a hot rarefied atmo- 
sphere (see Section 3.3). We let the initial torus evolve 
by solving non-radiative MHD equations, for the first 
elapsed time of t ~ 1 s, which corresponds to ~ 4.5 
times the Keplerian time at the center of the torus (at 
ro = 40i?s)- We then turn on the radiative terms in 
the basic equations and further solve the evolution of 
the torus for t ~ 10 s. Whereas the density normaliza- 
tion (which is the central density of the initial torus in 
the present study) can be taken arbitrarily in the non- 
radiative MHD simulations, simulations with the radia- 
tive terms do depend on the density normalization. That 
is, the density normalization controls the relative impor- 
tance of the radiative cooling. We will be able to repro- 
duce three distinct modes of accretion flow by changing 
the density normalizations (see Section 3.4). 

3.2. Code 

We numerically solve the set of RMHD equations us- 
ing an explicit-implicit finite difference scheme. The 
time step is restricted by the Courant— Friedrichs— Levi 
condition. The numerical procedure is divided into 
the following steps. In step I, the MHD terms 
are solved by the mo dified Lax— Wendroff scheme 
(jRubin fc BursteinI 119671 ). Equations (P)-© and ©, 
except the gas— radiation interaction terms of Equation 
(|3]), are solved in this step. The numerical p rocedure in 
this ste p is basically the same as that used in IKato et al.l 
(|2004al) , except that the equation of the internal energy 
of gas is solved in the present simulations, while they 
calculated the evolution of the total energy (internal and 
kinetic energies of the gas and magnetic energy). In step 
II, an integral formulation is used to generate a con- 
servative differencing scheme for the advection term of 
Equation Steps I and II are solved with the explicit 
method, while steps HI and IV are solved with the im- 
plicit method. The radiation energy and gas energy are 
updated simultaneously via the gas-radiation interaction 
in step HI. We consider the third and final terms of the 
right-hand side of Equation ([3]) as well as the terms on 
the right-hand side of Equation dH except for the radia- 
tive flux term in this step. In the final step (step IV), 
we update the radiation energy density via the radial 
radiative flux (the first term of the right-hand side of 
Equation (j4])). The radiation energy density is advanced 
again by the vertical radiative flux. In this step, the 
Thomas method is used for matrix conversion. 

Our method of calcul ation is an extensio n of those of 
MHD simulations fe.g.. IKato et al.l l2004aD. The MHD 
code was tested by the problems of the propagation of 
a two-dimensional MHD wave and of a one-dimensional 
shock, and applied to the simulations of the magnetized 
disks and jets (see also iKato et aTl l2004b[ ). We also 
performed the tests of gas— radiation interaction, one- 
dimensional and two-dimensional radiation front propa- 
gation, and sub- and super-critical shock, and radiation- 
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dominated shock, fin ding that the results a re in good 
agreement with those iTurner fc Stond ()2QQlf ) . 

3.3. Initial Conditions 

We initiahy set a rotating torus in hydrostatic balance. 
We employ a poly tropic equation of state, p = p^^^ with 
n = 3 and a power-law-specific angular-momentum dis- 
tribution, / = lo{r/ro)'' with Iq = {GMr^y/^/{r - Rs), 
ro = 40i?s, and a = 0.46. The initial density and gas 
pressure distributions are expressed as 



pt{r,z) = po 



1 



and 



Pt = po7eoV^^(^o,0) 



Pt{r,z) 



po 



l+l/n 



(18) 



(19) 



where po is the initial density at the center of torus, 
the parameter, eo, is set to be 1.45 x 10~^, and ?/^eff 
is the effective potential given by '0eff(^, ^) = V^l^,^) + 
0.5(//r)V(l-a). 

The initial magnetic fields in the torus are purely 
poloidal {B^ = 0). Their distribution is described in 
terms of the azimuthal component of the vector poten- 
tial, which is assumed to be proportional to the den- 
sity, A^p (X Pt. Other components are set to be zero, 
Ar = Az = 0. We initially set the plasma- /3, the ratio of 
gas pressure to magnetic pressure, to be 100. 

This initial torus is embedded in a nonrotating, hot, 
and rarefied corona with no magnetic fields. The density 
and pressure distributions of the corona are 



pc{r,z) = pi 



^{r,z)^GM/Rs 



CcGM/Rs 



and 



Pc{r,z) = pc{r,z) 



CcGM 



(20) 



(21) 



with pi = 10~^po and Cc = 1.0. This corona is initially 
in hydrostatic equilibrium. Our initial condition s are the 
same as those of model B in Kato et ajj (|2004bl ). except 
the density of the corona is smaller by a factor of 20. 

3.4. Models and Grids 

We need to assign the density normalization when 
starting RMHD simulations. In total, we calculate three 
models in the present study with po = Ig cm~^ (model 
A), lO^'^gcm"^ (model B), and 10"^gcm"^ (model C). 
We will show in Section 4 that our three models with 
high, moderate, and low density normalizations corre- 
spond to the slim disk, the standard disk, and the RIAF 
models, respectively. Note that we set the lower limit 
of the density to be p/po = 10~^^ for all models. Such 
a lower limit sometimes appears in the upper regions 
{z > 60Rs) around the rotation axis in model A, since 
the outfiowing velocity is so large. In models B and C, 
the density is rarely below the limit. 

Since we assume axisymmetry around the rotation 
axis, as well as refiection symmetry with respect to the 
equatorial plane, the computational domain can be re- 
stricted to one quadrant of the meridional plane. For 
models A and C, the number of grid points is (A/"^, Nz) = 



(512,512). The grid spacing which is uniform (Ar = 
Az = 0.2Rs)^ extends from 2Rs to 105i?s in the radial 
direction and from to 103i?s in the vertical direction. 
In model B, on the other hand, we use smaller spacing 
grids, Ar = Az = O.li^s, since the scale height of the 
disk is smaller (see below). The number of grid points is 
(Nr^Nz) = (1024,512), and the computational domain 
is 2Rs<r < 105Rs and < z < 51Rs. 

3.5. Boundary Conditions 

For the matter and magnetic fields, we adopt free 
boundary conditions at the inner and outer boundaries 
(r = 2Rs and 103i?s) and upper boundary {z = 103i?s 
for models A and C, z = 51 Rs for model B). That is, the 
matter can freely go out but not enter and the magnetic 
fields do not change across the inner, outer, and upper 
boundaries. If the radial component of the velocity is 
negative (or positive) at the outer (or inner) boundary, 
it is automatically set to be zero. The vertical component 
of the velocity is also set to be zero, when it is negative 
at the upper boundary. With respect to the disk plane 
{z = 0), we assume that p, p, Vr^ v^^ and Bz^ and are 
symmetric, while Vz, Br, and B^ are antisymmetric. 

The radiative fiuxes are assumed without using a gra- 
dient of the radiation energy density at the inner, outer, 
and upper boundaries. The vertical and radial compo- 
nents of the radiative fiuxes are set to be cEq at the outer 
and upper boundaries. At the inner boundary (r = 2Rs), 
we set the radial component of the radiative fiux to be 
zero except at < z < 2Rs. Since we also set the radial 
component of the advection of the radiation energy to be 
zero, the radiation energy neither leaves nor enters. At 
the boundary at r = 2Rs with < z < 2Rs, we assume 
Fq = —cEq, meaning that the radiation is swallowed by 
the black hole. In addition, we have Fq = at z = 0, 
since we impose a symmetric boundary condition at the 
equatorial plane. 

3.6. Prelusive Calculation 

As we have already mentioned in Section 3.1, we evolve 
the initial torus by solving non-radiative MHD equations 
for 1 s. Then, we solve only the MHD terms with grav- 
ity using step I in Section 13.21 In Figure [H we show 
the density distribution of the initial torus (top panel). 
Although the initial torus is in hydrostatic balance (see 
Section [33|) . the matter falls toward the black hole since 
the angular momentum is transported by the magnetic 
torque. Moreover, the magnetic energy dissipates and 
the flow is heated up. Thus, the initial torus also ex- 
pands in the vertical direction. The resulting density 
distribution at t = 1 s is shown in the bottom panel of 
Figure 1. Here, 512 x 512 grids (the grid spacing is 0.2Rs) 
are used. We have confirmed that the resulting density 
distribution does not change so much, even if we employ 
smaller mesh spacings, Ar = Az — O.li^s (1024 x 512 
grids). Assigning the density normalization, and setting 
the radiation temperature to be 10^ K in the whole re- 
gion, we start the RMHD simulations from the resulting 
structure given by the prelusive MHD calculations and 
go on performing them until t ~ 10 s. 

3.7. Mass Accretion Rates, Outflow Rates, and 
Luminosities 
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Fig. 1. — Color contour of the initial matter density distribution 
(top) and that obtained after 1 s non-radiative MHD simulations 
(bottom). 



The mass accretion rate is calculated at the inner 
boundary at r = 2Rs by 

r2Rs 

Mace = -2 / 27rrpvrdz, (22) 
Jo 

which is the mass passing through the inner boundary 
per unit time. The photon luminosity is calculated by 

Lph = 2 / 27rr {F§ + v,E^) dr, (23) 

at the height oi z = z^. The values of (rci,^c) are care- 
fully chosen so as not to include the contribution from the 
initial torus. We, hence, set (rci,2;c) = (25i?s7 60i?s) for 
model A, (30i?s, 30i?s) for model B, and (60i?s, 60i?s) 
for model C. In model C we employ relatively larger 
rci(= 60i?s), however, the contribution to the total emis- 
sion from the outer region is very small, since the density 
is very small there. The mass outflow rate and the kinetic 



luminosity are also calculated at z = as 

Mout = 2 / 2T:rpv,H {vr - v^sc) dr, (24) 

J2Rs 

Lkin = 2 ^ 27rr Qp^l^ "^zH {vr - Vesc) dr, (25) 

where H is the Heaviside step function; i.e., H{x) = 1 
for X > and H{x) = for x < 0, and vr[= Vr{r/R) + 
Vz{z/R)] is the i?-component of the velocity. We em- 
ploy (rc2,^c) = (60i?s, 60i?s) for models A and C, and 
(^c2,^c) = (30i?s, 30i?s) for model B. That is, what are 
meant by Mout and Lkin in the present study indicate the 
mass and kinetic energy transported upward only by the 
high- velocity outflow {vr > Vqsc) per unit time. In addi- 
tion, we plot the luminosity of trapped radiation, which 
is evaluated by 

r^Rs 

Ltrap = -2 / 27rr {F^ + VrEo) dz, (26) 
Jo 

with r = 2Rs. It implies the radiation energy swallowed 
by the black hole per unit time. 

4. RESULTS 

We first overview the simulation results in Section 14.11 
and then provide more detailed information in Sections 
14.21 - 14.41 for models A, B, and C, respectively. 

4.1. Overview of Simulated Flows 

The overall flow structures obtained by our RMHD 
simulations can be summarized in Figures [2]-[4j We first 
show perspective views of simulated flows in models A— C 
in Figure O Here, the color contours indicate the distri- 
butions of normalized density, p/po, time- averaged over 
6— 7 s for models A and C and over 9—10 s for model B. 
We find that a geometrically thick disk forms in models A 
and C, while a geometrically thin disk forms in model B. 
The streamlines indicated by the thick lines are overlaid 
in this figure. We find in all models that the gas near the 
equatorial plane is on a quasi-circular orbit around the 
central black hole, whereas the gas away from the equa- 
torial plane shows helical and outflowing motion. The 
helical motion means that the outflow material processes 
a substantial amount of angular momenta. 

The different dynamical properties of accretion flows 
in three models are more quantitatively shown in Figure 
[3l In the top panel of the figure, we plot the normalized 
mass accretion rate, Macc/(^E/c^), mass outflow rate, 
Mout/(^E/c^), photon luminosity, Lph/i^E, kinetic lumi- 
nosity, Lkin/^E, and trapping luminosity, Ltrap/^E- We 

also plot Lkin/^ph, ^trap/^ph, and Mout /Mace in the bot- 
tom panel. They are time- averaged over t = 5—7.5 s for 
models A and C and over 7.5 — 10 s for model B. For the 
calculation methods of these quantities, see Section [3771 
We find in model A that the photon luminosity exceeds 
the Eddington luminosity and that the trapping lumi- 
nosity is substantial, implying that the model A flow is 
supercritical flow. The thin disk in model B corresponds 
to the standard-type disk, since low scaleheight is a re- 
sult of efficient radiative cooling. By contrast, the pho- 
ton luminosity is negligible in model C, indicating that 
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Fig. 2. — Perspective view of inflow and outflow patterns near the black hole for models A, B, and C, from left to right, respectively. 
Also plotted are normalized density distributions (color) and streamlines, which are time- averaged over 6— 7 s for models A and C and over 
9-10 s for model B. 
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Fig. 3. — Top panel: mass accretion rate. Mace, and mass outflow rate, Mout, normalized by the critical rate, L-^/c^, 
kinetic, and trapping luminosities normalized by the Eddington luminosity, Lg. Bottom panel: the ratios of Mout/^aa 
^trap/^ph- All values are time-averaged over 5—7.5 s (models A and C) and 7.5—10 s (model B). 



and the photon, 
^kin/^ph, and 



the model C flow is RIAF. We also see a large value of 
^kin/^ph (> 1) only in model C. 

In Figure IH we show the radial profiles of the den- 
sity (top), the gas and radiation temperatures (middle), 
and the radial velocities (bottom) around the equato- 
rial plane. They are time-averaged over t = 7.5—10 s. 
As expected, the gas temperature is highest in model C 
(RIAF), while it is lowest in model B except at radii less 
than ~ 5i?s- The decoupling of the gas and radiation 
temperatures occur in model C entirely and in model B, 



inside ~ 7Rs. 

To summarize, we could reproduce three distinct states 
of accretion flow with the same code but by changing the 
density normalization. In the subsequent subsections we 
give more information individually for models A— C. 

4.2. Model A 
4.2.1. Overview 

Model A, with a high density normalization {po = 
l.Ogcm"^), corresponds to the two-dimensional RMHD 
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Fig. 4. — Radial profiles of the density (top), the gas and radiation temperatures (middle), and the infiow velocity normalized by the 
local Keplerian velocity (bottom) near the equatorial plane for models A (z = O.IRq, black), B (z = 0.05Rs, red), and C (z = O.IRq, blue). 
All values are time-averaged over 7.5—10 s. 



version of the slim dis k model (jAbramowicz et aTl 
[1988'; "Watarai et al.' "2000) but with significant outflow 
(JTakeuchi et al. 2009). The geometrically thick disk is 
supported by the radiation pressure (see Section 14.2.21 
and Section l4.2.3p . We also flnd in Figure H] that the 
disk consists of very dense and moderately hot (10^~^ 
K) matter. The gas temperature is comparable to the 
radiation temperature due to the effective gas— radiation 
interaction. The inflow velocity {—Vr) tends to increase 
inward only in the vicinity of the black hole, r<7Rs^ 

and it is very small or even negative (i.e., outflow) around 
r lORs (see the bottom panel in FigureH]). Basic prop- 
erties of model A are roughly consistent with those of the 
slim disk. However, the slopes of the radial proflles of the 
density and temperatures do largely differ from those of 
the slim disk. (Note that such deviations from the con- 
ventional disk models are also found in models B and C.) 
These discrepancies might be due to the limited calcula- 
tion time, the restricted computational domain, and the 
influence of the initial conditions. Detailed comparison 
is left for future work. 

In Figure [3l we flnd that the mass accretion rate 
onto the black hole is much larger than the critical 
value. Mace ^ 10^Le/c^. The photon luminosity is 
^ph 1.7 Le, which is calculated based on the vertical 
component of the radiative flux within the polar angle 
of l9 = tan-i(25i?s/60i?s) ^ 23°. (Note that we as- 
signed (rci,^c) = (25i?S7 60i?s) in Equation (|23]) .) Since 
(^ci,^c) = (25i?s, 60i?s) lies o n the line of the photo- 
sphere, z/r ~ 2.4 (see Section l4.2.2p . the photon lumi- 
nosity calculated with these values is equal to the sum of 
the radiation energy released at the photosphere of the 
inner part of the disk, r < 25Rs^ per unit time. Also, 
since rd = 25Rs is smaller than the position of the ini- 



tial torus, r = 40i?s, we can avoid possible influence by 
the initial conditions. Here we note that the calculated 
luminosity (i^ph) may be underestimated, since the con- 
tribution from the outer disk (r >25i?s) was ignored. For 
better calculations of Lph, we need simulations with a 
larger computational domain, in which we set the initial 
torus at a much larger radius. 

In addition, the photons mainly escape along the ro- 
tation axis, producing the mildly collimated emission 
(note that the vertical component of the radiative flux 
is much larger than the radial one around the rotation 
axis, z/r > 2.4). Thus, the apparent luminosity of the 
flow would exhibit the strong viewing-angle dependence. 
We will discuss this issue in Section O 

Figure [3] also indicates that the trapping luminosity 
is slightly larger than the Eddington luminosity (top 
panel), and is comparable to the photon luminosity (bot- 
tom panel). This means that the large number of pho- 
tons generated inside the disk is swallowed by the black 
hole together with accreting matter, without many being 
released from the disk surface. 

Although we do not represent the rotation velocity in 
Figure IH the matter approximately rotates with the Ke- 
plerian velocity (within 10%) near the equatorial plane. 
The outflow is driven by the radiation force (see Sec- 
tions 14.2.21 and I4.2.3P . The mass-outflow rate is a few % 
of the mass accretion rate (see Figure [3]). We also flnd 
^kin <0.1Lph in this flgure. This implies that the su- 
percritical flows in model A release energy via radiation 
rather than via outflows. 

4.2.2. Two-dimensional Structure 

Figure[5]illustrates various aspects of the flow structure 
for model A. The contours of the density, the gas temper- 
ature, and the plasma- /3 are shown in panels (a), (b), and 
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Fig. 5. — Two-dimensional distribution of the various quantities for Model A: (a) the density overlaid with the velocity vectors, (b) the gas 
temperature, (c) the plasma-/3, (d) the magnetic energies via the toroidal component of field, (e) the same but of the poloidal component, 
(f) the magnetic pitch, (g) the radiation energy, (h) the ratio of the radiation energy to the sum of the gas and magnetic energies, (i) and 
the ratio of the gas temperature to the radiation temperature. All values are time-averaged over t = 6—7 s. The white and black arrows 
in panel (a) indicate the velocity vectors whose magnitude exceed the escape velocity. The dashed line in panel (b) is the photosphere, at 
which the optical thickness measured from the upper boundary is unity. The arrow in panel (g) shows the radiative flux vector. 
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(c), respectively. Here, the velocity vectors are overlaid 
in panel (a) with arrows, in which the white arrows indi- 
cate that the velocity exceeds the escape velocity defined 

1/2 

as Vesc = {2GM/R) ' . The magnetic energies of the 
toroidal and poloidal components are separately plotted 
in panels (d) and (e), respectively, where the poloidal 

1 /2 

component is given by Bp = (5^ + • The mag- 

netic pitch, \B^p\/Bp, is shown in panel (f). Panel (g) 
indicates the radiation energy, and the ratio of the radi- 
ation energy to the sum of the gas and magnetic energies 
is plotted in panel (h). Panel (i) shows the ratio of the 
gas temperature to the radiation temperature, where the 
radiation temperature is calculated as Trad = (^o/<^r)^^^- 
All values are time- averaged over t = 6—7 s, and both 
axes are normalized by the Schwarzschild radius. 

We find in panel (a) that the flow is divided into two re- 
gions: the disk region around the equatorial plane (char- 
acterized by white and red colors) and the outflow region 
above the disk region (blue region). The boundary be- 
tween the two regions is approximately along the line of 
z/r ^ 2. The matter slowly accretes onto the black hole 
through the disk region with velocity much less than the 
free-fall velocity (see small black arrows). This panel also 
shows that the outflow velocity highly exceeds the escape 
velocity (see white arrows). The photosphere, at which 
the optical thickness measured from the upper calcula- 
tion boundary is about unity, roughly corresponds to the 
green region, i.e., z/r ^ 2.4. 

We find that the hot, rarefied corona appears above the 
supercritical disk. Panel (b) indicates that the outflowing 
matter is very hot (^10^ K), whereas the disk consists 

of relatively hot gas (~ 10^~^ K). The gas density is 
much smaller in the outflow region than in the disk region 
(see panel (a)). Such a hot outflowing cor ona above the 
supercr itical disk has also been reported bv lOhsuga et al.l 
(|2005b), and is though to affect the spectra due to the 
Comptonization (Kawashima et al. 2009). 

Panel (c) shows that the outflow is surrounded by the 
low plasma- /3 region (blue), in which we find plasma- /3 
< 1. In such a low plasma- /3 region, the toroidal compo- 
nent of the magnetic fields, which seems to be amplified 
via the differential rotation of gas, is dominant over the 
poloidal one (see panels (d) and (e)). This feature is 
also clearly demonstrated in panel (f), in which we see 
\By:^\/Bp > 10 around the outflow region. That is, the 
outflow is surrounded by the regions with strong toroidal 
magnetic fields. Inside the outflow, on the other hand, 
the vertical component of the magnetic fields is dominant 
over the other components. 

Such a magnetic stru cture is quite remi n iscent 
of magnetic-tower jets (iLynden-Bell fc Boilvl 119941 : 
iLvnden-Belil 119961 : iKato et al.l l2004b[ ). However, there 
is one big difference; while the magnetic tower jets are 
accelerated via the magnetic-pressure force, the outflow 
in our model A is powered by the radiation force. This 
radiatively driven outflow is collimated by the Lorentz 
force. Thus, on the basis of global RMHD simulations 
we proposed a novel jet model; the radiation-pressure 
driven and magneti cally collimated outflow (see also 
iTakeuchi et"alll2010f ). 

Note that there must be something to prevent the 
field from expanding sideways, thereby realizing the high 




Fig. 6. — Vertical structure at r = lOi^s for model A. Top panel: 
the density (red) and vertical component of the velocity normal- 
ized by the escape velocity (black). Middle panel: the gas and 
radiation temperatures. Bottom panel: the vertical components of 
the gravity (dotted line), the radiation force (green filled circle), 
and the magnetic-pressure force (blue filled circle). The sum of the 
vertical forces by the radiation, the magnetic pressure, the mag- 
netic tension, and the gas pressure is plotted by the black solid 
line (total). Here, the gas-pressure force (red filled circle) and the 
magnetic-tension force (blue open circle) are not represented, since 
they are so small or negative. All values are time-averaged over 
6-7 s. 

magnetic energy density in the cylindrical region around 
the black hole. In our simulations it is the radiation- 
pressure force by geometrically thick flow surrounding 
the magnetic tower that confines the magnetic tower. 
The magnetic field can thus naturally evolve to form 
such a high concentration from its initial configuration 
(Take uchi et al. 2010). 

In panel (g), it is found that the radiation energy is 
enhanced in the disk region. The disk is not only ge- 
ometrically thick but also optically thick. Because of 
numerous scatterings, the photons generated in the disk 
cannot easily escape from the disk surface. As a conse- 
quence, a large number of photons accumulate in the disk 
region, leading to the high radiation energy density. The 
radiation energy is predominant over the gas energy, and 
also over the magnetic energy in the whole region (see 
panel (h)). The ratio of the radiation energy to the sum 
of the gas and magnetic energies is >10^ in the outflow 
region and ^10 in the disk region. 

Panel (i) indicates that the gas temperature is nearly 
equal to the radiation temperature within the disk (blue). 
The absorption opacity is large in the disk region, since 
it is sensitive to the density (ex p^). Thus, the effective 
gas— radiation interaction (emission/absorption) leads to 
^gas ^ Trad' In coutrast, the gas temperature is much 
higher than the radiation temperature in the outflow re- 
gion, where the radiative cooling is not efficient because 
of the small absorption opacity (small density). Thus, 
the gas temperature is kept high (see the next subsec- 
tion for more quantitative descriptions). 

4.2.3. Vertical Structure 

In the top panel of Figure [6l we plot the vertical pro- 
files of the density (red) and the vertical component of 
the velocity (black) for r/Rs = 10. Here, they are time- 
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Fig. 7. — Vertical profiles of the energy dissipation rate (black) 
the gas (red), radiation (green), and magnetic (blue) energies at 
r = 10i?s- They are time- aver aged over 5—7.5 s. The radiation 
energy in model C is not plotted, since it is too small. 

averaged over t = 5—7.5 s. We find that the density 
decreases with an increase of z. The upward velocity in- 
creases as z increases and it exceeds the escape velocity 
at z>20Rs. The sonic point is around z = 20i?s- While 
the gas temperature is comparable to the radiation tem- 
perature due to the effective gas— radiation interaction 
at z<15Rs^ we find Tgas ^ T^ad at the upper region of 
z>20Rs, since radiative cooling is inefficient because of 

smaller density (emissivity) (see the middle panel). 

As we have mentioned in Section |4.2.2[ the supercrit- 
ical disk ejects high- velocity outflows driven by the ra- 
diation force. This is clearly shown in the bottom panel 
of Figure [6l in which we plot vertical components of the 
gravity (dotted line), the radiation force (green circle), 
the magnetic-pressure force (blue circle). The solid black 
line (total) in the bottom panel indicates the sum of the 
vertical forces by the radiation, the magnetic pressure, 
the magnetic tension, and the gas pressure. (The gas- 
pressure force and the magnetic-tension force are so small 
that they do not appear in this figure.) While the radia- 
tion force is nearly balanced with the gravity at z<20Rs 

(disk region), it largely exceeds the gravity in the upper 
regions at z>20Rs, producing the high- velocity outflows. 
The radiative flux is roughly estimated to be cEq in the 
optically thin, outflow region above the disk, whereas 
it is reduced to be cEq/t in the optically thick disk re- 
gion. Hence, the radiation force suddenly increases above 
around the disk surface, z ~ 20Rs^ although the radia- 
tion energy density itself (Eq) is smaller in the outflow 
region than in the disk region (see panel (g) of Figure [5j). 
In the case of the supercritical flows, the radiation force 
supports the disks and accelerates the outflows. The 
forces via the magnetic pressure, the magnetic tension, 
and the gas pressure is too small to influence the accel- 
eration of the outflow. However, the magnetic pressure 
and the magnetic tension (Lorentz force) works in a di- 
rection para llel to the disk p lane and thus collimates the 
outflow fsee iTakeuchi et al.ll2QlQ[ ). 

4.2.4. Dissipation Rate 



In Figure El we plot the vertical profiles of the dissipa- 
tion rate, AtttjJ'^/c^ (black), as well as those of the gas 
(red), radiation (green), and magnetic energies (blue) for 
model A (see upper four lines). They are time- averaged 
over 5—7.5 sec. All the values tend to decrease with an 
increase of z. In the region of z<20Rs^ we find that the 
gradient of the dissipation rate is smaller than that of the 
gas energy and that of the radiation energy. It implies 
that the traditional a- viscosity model does not precisely 
describe the vertical structure of the disk, since it states 
that the dissipation rate is proportional to the gas (or ra- 
diation) pressure. We find that the profile of the dissipa- 
tion rate is rou ghly on parallel with that of the m agnetic 
energy (see also lTurner et al.|[2QQ3l : lHirose et al.| [2QQ6. for 
the cases of local RMHD simulations). 

The radiation energy shows a rather flat distribution 
above z>20Rs. Since the radiation energy is transported 
via the diffusion in the disk region, while the photons 
freely go out above the disk, the slope of the profile of 
the radiation energy changes across the disk surface at 
z - 20Rs. 

In figure [71 we also find that, while the dissipation rate 
is enhanced near the equatorial plane in both models, 
the magnetic energy suddenly decreases toward z = 0. 
However cautions should be taken here, since this might 
be caused by the particular boundary conditions with re- 
spect to the equatorial plane, where we require that Br 
and are antisymmetric. Such a condition enhances 
the reconnection, leading to the decrease of the mag- 
netic ener gy and t he increase of the dissipation rate. In 
fact, Hiros eet al.l (2006) showed that the magnetic en- 
ergy and the dissipation rate are nearly constant across 
the equatorial plane by simulations without employing 
the equatorial-plane symmetries. 

Although the traditional a-viscosity model does not 
precisely describe the vertical dissipation profile, the r—(f 
component of the shear-stress tensor at around the equa- 
torial plane is roughly proportional to the pressure when 
we examine the time va riations of these two quantities. 
In the previous paper (jOhsuga et al.l l2QQ9f ) , we inves- 
tigated the time evolution of the magnetic torque of 
z<3Rs^ and demonstrated that the torque is roughly 
proportional to the total pressure but with some scat- 
ters. 

4.2.5. Gaussian or Polytropic? 

Finally, we compare the simulated vertical structure 
with that obtained by simple analytic models on the as- 
sumption of the hydrostatic balance. If the disks are ver- 
tically isothermal, the energy (pressure) profile is given 
by a Gaussian profile of exp (— 2;^/2i^^) with H being 
the disk half- thickness. Such a function with H = 2.5Rs 
roughly reproduces the energy distributions of the flow 
in model A, as is seen in the top panel of Figure [HI which 
represents gas and radiation energies of model A normal- 
ized by the values at z = 0. If we use a polytropic rela- 
tion {p oc with N being the polytropic index) 
instead of the isothermal assumption, on the other hand, 
the vertical hydrostatic balance leads that the energy is 

proportional to a function of (l — z'^ /2H'^^^^^ . Such a 
function with H = 8i?s can also give good fits to the en- 
ergy profiles of model A (see the upper panel in Figure 
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Fig. 8. — Gas and radiation energies of model A (top) and gas 
energy of model C (bottom) normalized by the each value at z = 
0. Solid and dotted lines indicate the Gaussian and polytropic 
relations, respectively. 

[8]). Here, we employ = 3 for model A since the disk 
is radiation-pressure-dominated. We note that the disk 
half-thickness at r = lOi^s is ~ (cs/'^Kep) ^ 5i?s and 
are in rough coincidence with the above values, where cs 
is the sound velocity at (r, = (10i?s,0). 

4.3. Model B 

4.3.1. Overview 

As is already mentioned, accretion flows at r>7Rs in 

our model B (po = 10~^gcm~^) resemble the standard 
disk. The efficient radiative cooling leads to the forma- 
tion of a geometrically thin but optically thick disk (see 
Figure [2]). The disk in model B is moderately dense as 
shown in the top panel of Figure IH The gas temperature 
of ~ lO^K at r>7Rs is close to the value predicted by the 

standard disk model (see the middle panel of Figure |4]). 
However, the disk is truncated around r ~ 7Rs (which 
we call the truncation radius). Because of small density 
of the accretion flow near the black hole, the emissivity 
is greatly reduced, which is responsible for a decoupling 
of the gas and radiation temperatures, Tgas ^ Trad at 
r <5i?s- The bottom panel shows that the inflow veloc- 
ity {—Vr) increases as the flow approaches the black hole. 
The rotation velocity is very close to the Keplerian ve- 
locity (within 10%), although we do not represent it in 
Figure HI 

We find in Figure [3] that the photon luminosity is well 
below the Eddington luminosity, ~ 2 x 10~^Le, for the 
mass accretion rate is ~ 5 x 10~^ L^/c^. In this model, 
the photon trapping effect is negligible. The energy con- 
version efficiency, rj = Lph/MaccC^ ^ 0.04, is the largest 
in three models. This is also consistent with the pre- 
diction of one-dimensional accretion disk study, whereby 
the radiative cooling is more effective in the standard 
disk model than the slim disk model and RIAF. How- 
ever, the simulated value of r] is smaller than the pre- 
diction of the standard disk model, ~ 0.1. This result 
might be caused by the disk truncation around r ~ 7Rs^ 
within which the density is too small for emission to be 
efficient (will be discussed in Section [5^ . If we were to 



employ higher density normalization, both the photon 
luminosity and the conversion efficiency might increase, 
since then the truncation radius would decrease or even 
disappear. Such a disk would be more consistent with 
the standard disk model. 

Whereas the disk wind is not predicted by the standard 
disk model, a significant amount of matter is blown away 
in our model B. As we have already mentioned in Section 
14.11 the matter above the disk rotates and goes upward. 
In other words, outflow material possesses substantial 
angular momentum, taking it away from the underly- 
ing disk. Figure [3] illustrates that time-averaged mass- 
outflow rate is about 10~^ Le/c^ (see the top panel), 
which is a few % of the mass accretion rate (see the bot- 
tom panel). We also find in the bottom panel that Lkin 
is a few % of Lph- As is the case in model A, the disk 
in model B releases energy via radiation rather than via 
outflows. Here we note a possibility that the outflowing 
matter with high velocity might originate not from the 
outer cold part of the disk {r>7Rs) but from the inner 

hot part of the disk {r<bRs). The matter ejected from 
the outer cold part of the disk might produce low- velocity 
outflow with a wide opening angle. At this moment it is 
technically very difficult to identify the launching points 
of the high- and low- velocity outflows. Note that there 
is no clear evidence for the collimation of the outflows in 
model B. 

4.3.2. Two-dimensional Structure 

Two-dimensional structure of the flow in model B is 
presented in Figure [9l Each physical quantity is time- 
averaged over t = 9—10 s. Panel (a) clearly shows that a 
geometrically thin disk is located at around the equato- 
rial plane (indicated by the white and red colors). The 
disk thickness is very small and is everywhere less than 
^ IRs. The gas of the disk effectively cools by emitting 
radiation, as shown in panel (b), except at r<bRs. The 

gas temperature is very high above and below the disk, 
since the low density leads to a low cooling rate (emis- 
sivity). Such high-temperature and low-density matter 
forms a corona surrounding the cold disk and Compton 
upscatters seed photons generated within scale height the 
disk. 

Panel (c) shows that a relatively low plasma- /3 region 
appears along the line of |z|/r ~ 2 (green and light blue). 
We find that the matter goes upward in this region (see 
vectors in the blue region in panel (a)). Although the 
time-averaged velocities are less than the escape veloc- 
ity, the upward velocity intermittently exceeds the escape 
velocity. As a result, the time- averaged mass-outflow 
rate is Mout IQ-'^Le/c^ (see Figure [S]). Panel (f) 
shows that this outflow is surrounded by the regions of 
\By:^\/Bp > 1, while the poloidal (vertical) component of 
the magnetic fields is dominant at the very vicinity of the 
rotation axis. These features can be seen in panels (d) 
and (e). The toroidal field is enhanced mainly near the 
disk, |z| <10i^s (see panel (d)). The poloidal component, 

in contrast, tends to be enhanced near the black hole (see 
panel (e)). Again, something like magnetic tower struc- 
ture appears in model B, as well. Here we note that the 
higher plasma- /3 and the larger magnetic pitch at the 
outer region {r>20Rs and \z\<15Rs) are thought to be 
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Fig. 9. — Same as Figure [5] but for model B. Each quantity is time- aver aged over t = 9—10 s. Note the different plot area from that of 
Figure [5] 
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Fig. 10. — Same as Figure [6] but for model B. Each quan- 
tity is time- aver aged over t = 9—10 s. Here, the radiation force 
(green filled circle) and the magnetic-tension force (blue open cir- 
cle) are not represented, since they are so small or negative. Here 
we note that the disk region contains 5—6 grid points (excluding 
the point at the boundary). Although this number is not big, it 
can marginally resolve the vertical structure. 

under the influence of the initial torus. 

We find in panel (g) that the radiation energy is en- 
hanced within the disk (white and red), since the disk is 
optically thick except at r<7Rs. The optical thickness of 
the disk atmosphere is, by definition, smaller than unity, 
implying that the photons emitted at the disk surface 
can freely go out. Thus, the radiation energy is small 
above and below the disk (blue). In contrast with model 
A, the radiation energy is smaller than the gas energy 
and the magnetic energy (see panel (h)). 

Panel (i) shows Tgas ~ T^ad in the disk at r>7Rs 

(blue), and Tgas ^ ^rad in the other regions (white and 
red). This is because the gas— radiation interaction is ef- 
fective only in the dense region as we have already men- 
tioned. Such a feature is similar to that in model A, 
though the disk is geometrically thin in model B. 

4.3.3. Vertical Structure 

Figure [10] is the same as Figure [6] but for model B 
and t = 7.5—10 s. In the bottom panel, the red circles 
indicate the gas-pressure force. In the top panel, the 
high-density regions of z < Rs correspond to the thin 
disk as we have shown in Figure [9l The density slowly 
decreases with an increase of z above the disk. This is 
because the gas temperature is very high and so is the 
scale height in the low-density region dit z > Rs (middle 
panel). We find Tgas ^ ^rad in the disk region. 

We find that sum of the gas and magnetic-pressure 
forces roughly balance with the gravity at z>3Rs (see 

the bottom panel). The time-averaged upward velocity 
is much smaller than the escape velocity in this region 
(see the top panel). Note, however, that the high- velocity 
outflows occasionally appear. 

How is the matter ejected from the disk? It is mainly 
by the magnetic-pressure force. The bottom panel shows 
that the hydrostatic balance breaks down around the disk 
surface, z = 1 — 3Rs^ where the magnetic-pressure force 
exceeds the gravity, leading to the mass ejection from 
the disk surface. In contrast with model A, the radiation 



force is negligible. 

Although we used small grid spacings, Ar = Az = 
O.li^s, in our simulations for model B, we note that the 
simulations of even higher resolution are required to in- 
vestigate the detailed structure of the geometrically thin 
disk. Such calculations will be performed in a future 
study. 

4.4. Model C 

4.4.1. Overview 

Flows in model C (po = 10~^gcm~^) correspond to 
the RIAF. The geometrically thick disk is supported by 
gas pressure i n coo pera tion w ith the magnetic pressure 
(see Sections 14.4.21 and I4.4.3p . As shown in Figure HI 
the flow is entirely occupied by very hot and rarefied 
plasma. Since the density is too low for radiative cooling 
to be effective, a decoupling of the gas and radiation 
temperatures, Tgas ^ ^rad, occurs in the whole region. 
The bottom panel in Figure |4] shows that the accretion 
velocity is relatively high at r<10Rs in comparison with 

those of the flows in models A and B 

Since the optical thickness is very small (r ~ 10""^), 
the flow in model C is radiatively inefficient and faint. 
The photon luminosity is much smaller than the Edding- 
ton luminosity, Tph <C 10~^Te (see Figure [3]). The en- 
ergy conversion efficiency in model C {rj ~ 10 ~^) is much 
smaller than those in models A and B. 

The helical streamlines seen in Figure [2] indicate the 
formation of jets around the rotation axis. The mass- 
outflow rate is about 10 % of the mass accretion rate 
(see Figure [3]). Interestingly, we find Tkin > lO^Tph in 
model C (cf. Tkin < O.lTph in models A and B). Hence, 
the radiatively inefficient disk (Mace <^ Le/c^) loses the 
energy more via jets rather than via radiation. 

4.4.2. Two-dimensional Structure 

Figure [TT] shows the two-dimensional structure of the 
flow in model C. This figure is the same as panels (a) — (f) 
of Figure [5] but for model C. The white and red colors 
in panel (a) indicate a geometrically thick disk (disk re- 
gion). The entire zone can be divided to the disk and 
outflow regions by the line of z/r ^ 2. The gas tempera- 
ture highly exceeds 10^ K in the whole region (panel (b)). 
Here we note that we might overestimate the gas temper- 
ature, since the cooling via the Compton scattering is not 
taken into consideration. Also, the electron temperature 
might deviate from the ion temperatures, although we 
treat the matter as a one-temperature plasma. The disk 
is mainly supported by the gas pressure. The gas energy 
is dominant over the magnetic pressure in the disk region 
(see panel (c)). 

In panel (a), we also find that the matter with rel- 
atively small density is blown away above the disk, i.e., 
z/r>2. (The outflow region is characterized by the green 

and yellow colors.) The time- averaged outward velocity 
exceeds the escape velocity only in the region of z>70Rs 

in panel (a) (white vectors). Note, however, that this fig- 
ure only shows the time average, while the simulations 
exhibit significant time variations in vr. In fact, we find 
that the high- velocity outflows {vr > Vqsc) occasionally 
appear at z>10Rs near the rotation axis. The outflow is 
mainly accelerated via the magnetic pressure, in cooper- 
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Fig. 11. — Same as panels (a) — (f) of Figure O but for model C. Here, note that there is no photosphere, since the flow is optically thin. 
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Fig. 12. — Same as Figure[6]but for model C. Here, the radiation 
force (green filled circle) is not shown, since they are negligibly 
small. Note that the gas pressure exceeds the total pressure at 
small height, z < 2Rs (see the bottom panel). This is because the 
vertical force by the magnetic pressure is negative there. 

ation with the gas pressure (we will discuss below). The 
radiation force is negligible. 
Panel (c) indicates plasma-/^ <1 in the outflow region, 

whereas /3 > 1 in the disk region. The toroidal mag- 
netic fields are mainly amplified in the disk region, while 
poloidal fields strengthen around the rotation axis (see 
panels (d) and (e)). Hence, the magnetic pitch, |5(^|/5p, 
tends to increase as an increase of r. Panel (f) shows that 
the pitch becomes large around the outer edge of the 
outflow region located long the line of z/r ~ 2 (see the 
yellow area). This implies that the magnetic field lines 
are strongly coiled around the spine of the outflow; that 
is, helical magnetic fields form. Thus, the outflow in our 
model C seems to be magnetic-tower jets, which have 
been re ported by the three-dime nsional MHD simula- 
tions bvlKato et a l. (2004b; see also lLvnden-Bell fc Boilvl 
11994"; Lynden-Beh 1996, for the original proposal). 

Since the disk is optically thin, the photons can freely 
escape from the disk, producing a quasi-spherical distri- 
bution of the radiation energy density, as is clearly shown 
in panel (g). Again, the radiation force is negligible due 
to the small radiation energy and small opacity. 

4.4.3. Vertical Structure 

Figure [12] is the same as Figure [6] but for model C. The 
density and the gas temperature monotonically decrease 
with an increase of z. The open blue circles indicate 
the vertical component of the magnetic-tension force. 
We find that the geometrically thick disk at z<20Rs is 
mainly supported by the gas pressure, achieving the hy- 
drostatic balance (see the bottom panel). The magnetic- 
pressure force is dominant over the gas-pressure force 
above the disk, z>20Rs. In this region, the total force 

(black line) exceeds the gravity, implying that the matter 
is accelerated upward. That is, the magnetic-pressure 
force in cooperation with the gas-pressure force drives 
the outflows above and below the disk. The outflow is 
thus accelerated upward; its velocity exceeds the sound 
velocity at z ~ 30Rs and reaches the escape velocity at 
z ^ bORs. Both the radiation force and the magnetic- 
tension force are found to be negligible. 



Similarly to model A (see Figure [6]), we find in model 
C that the disk is in hydrostatic balance and the mat- 
ter is blown away from the disk surface. However, the 
radiatively driven outflows in model A is more power- 
ful than the magnetically driven outflows in model C 
(The kinetic luminosity, Lkin, is larger in model A than 
in model C). While the radiation force largely exceeds 
the gravity in model A, the upward force only slightly 
exceeds the gravitational force in model C (see the bot- 
tom panel in Figure [T2]) . As a result, low-luminosity 
disks with Mace <C /c^ cannot produce powerful high- 
velocity outflows in a steady fashion, but only occasion- 
ally. 

4.4.4. Dissipation Rate 

The three lower lines in Figure [7] represent the verti- 
cal profiles of the dissipation rate (black), gas (red) and 
magnetic energies (blue) for model C, which are time- 
averaged over 5—7.5 s. Here, we note that the radiation 
energy is never appreciable in model C, so it is not plot- 
ted. We find that the slope of the profile of the dissipa- 
tion rate is nearly equal to that of the magnetic energy, 
and flatter than that of the gas energy. Such a feature 
implies that the traditional a-viscosity model does not 
precisely describe the vertical flow structure, as was al- 
ready noticed in model A. As we have discussed in Sec- 
tion I4.2.4[ the boundary condition with respect to the 
equatorial plane seems responsible for the enhanced dis- 
sipation rate, leading to the drop of the magnetic energy 
at z ~ 0. 

4.4.5. Gaussian or Polytropic? 

In Section I4.2.5[ we have mentioned that the energy 
distributions of the flow in model A agrees with the pro- 
files of the disks in hydrostatic balance (isothermal and 
polytropic). However, we cannot fit the profile of the 
gas energy of model C by the Gaussian profile (with 
H = 3Rs) nor by the polytropic relation (with H = 7Rs, 
see the bottom panel in Figure [8]). Here, = 1.5 is used 
because the gas pressure is predominant for the disk of 
model C. Even if we change fitting is not success- 
ful. The drop of the gas energy is more rapid at small 
z {<4:Rs) and steeper above. Such deviations may be 
caused by the particular magnetic energy profile which 
is not simply proportional to gas pressure. Finer-mesh 
calculations are necessary to confirm if this is really the 
case. 

5. DISCUSSIONS 

5.1. Comparison with Observations: Outflow 

Our simulations demonstrate that powerful outflows 
appear in the super- or near-critical flows (model A). 
This result seems to be consistent with the observa- 
tions of bright black hole objects with high Eddington 
ratio, such as narrow- line Seyfert 1 galaxies (NLSls) 
and microquasars. The NLSls, which are thought to 
contain relatively less massive black holes and thus to 
be high Lph/I/E systems, are usually radio-quiet. How- 
ever, som e of the m have been reported to be radio- loud. 
iDoi et aP ()2006l ). for example, concluded by very long 
baseline interferometry observations that a radio-loud 
NLS l, J094857 . 3+00 2225, possess relativistic jets (see 
also lZhou et al.ir2003f ). By contrast, there are plenty of 
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such sources in our Galaxy. Microquasars, such as GRS 
1915 +105, also sho w relativistic jets in high- luminosity 
state (jFender et al .1120 04) and SS433 is another good ex- 
ample for supercritical accretion flow with powerful jets. 

Our model A may explain the nebulae (or bubbles) 
around ultraluminous X-ray sources (ULXs) or micro- 
quasars. iPakull et al.l ()2010j i proposed by X-ray/Optical 
observations that the large nebula, S26, is produced by 
energy injection via powerful jets, of which the kinetic 
luminosity is estimated to exceed the photon luminos- 
ity, I/kin > -^nh: by three orders of magnitude (see also 
iPakuh fc Grisi[2008f ). The supercritical flows will be ob- 
served as low luminosity objects (Lph < i^kin) near the 
edge-on view, since the emission is mildly collimated and 
since the radiation from the innermost region can be ob- 
scured by the outer disk. If this is the case, the large 
^kin/^ph ratio of S26 can be explained by our model A. 
By solv ing!: the radiation transfer equation, iQh suga et al.l 
(|2005b[ ) demonstrated that the apparent luminosity in 
the edge-on view could be more than 10 times smaller 
than that in the face-on view. In this calculation, large 
luminosity comes from the middle part of the disk at 
r>100Rs by the limitation of the computational box. 
If we can make simulations with a much more extended 
computational box, an even smaller value of Lph is ex- 
pected at the edge-on view. Since Lkin is on the order 
of Le regardless of the viewing angle, the large Lkin/^ph 
is, hence, expected. 

Our model B shows the disk wind and the geometri- 
cally thin disk, which was not predicted in the frame- 
work of the standard disk model. The recent X-ray 
observations of black hole binaries (BHBs) reveal the 
presence of the blueshifted absorption lines, meaning 
the occurr ence of the d i sk win d from the geometrically 
thin disk. iMiller e t all (l2006allbl) concluded by the ob- 
servations of GRO J1655-40 and a black hole candi- 
date, H1743-322, that the matter with p - lO^^gcm"^ 
ejected from the geometrically thin disk with the speed 
of 500 - 1600 km s-^ The outflow velocity of the X- 
Ray transient, 4 U 1630-472, has be en also reported to 
be - 10^ kms-^ (jKubota et al.ll2007[ ). The present study 
can account for such observations. Our model B shows 
that the gas with p ^ 10~^ gcm~^ goes out at the speed 
of ^ lO^kms"^ (see Figures [9l and [T0|) . 

We find the magnetic-tower jet in model C. In contrast 
with models A and B in which we find Lkin < ^ph, the 
kinetic luminosity largely exceeds the photon luminosity 
in model C, Lkin ^ ^ph (see Figure [3]). That is, the 
radiatively inefficient flow releases the energy mainly via 
the jet. This result can be understood by the predic- 
tions, whereby the kinetic luminosity is proportional to 
the mass accretion rate, Lkin oc M^cc but L^oc M\ 



m tne radiatively memcient regime |-bender et al.l I^UU4J : 
iNaravan fc Yil[T995[ ll^o et al.ll2008L see also Figure [3]). 

5.2. Comparison with Observations: Spectra 

Our model A might resolve the outstanding issue 
regarding the central engine of ULXs. As for their 
central engine, two competing models have been dis- 
cussed over the past decade: subcritical accretion 
onto intermediate-mass black holes wi th the black hole 
mass greatly exceeding lOOM© (e.g., Ma kishima et al.l 
I2OOO; Miller et al...2004) and supercritical accretion onto 



stellar- mas s black holes with mass of 3 — 2OM0 (e.g., 
IKing et al.l 120011 : IWatarai et all I2001D . In addition, 
there is the possibility that some ULXs may be pow- 
ered by stellar (but not stellar-mass) black holes in the 
30 — 8OM0 range, formed in low metallicity environ- 
ments, with near-critical or slightly supe r critical accre- 
tion flows (e.g., IZam pieri fc Robert si 120091 : iMapelli et al.l 
I2OIOI : iBelczynski et al..,.2010f ). It is thus crucially impor- 
tant to show theoretically how much (apparent) lumi- 
nosities can be achieved by supercritical accretion. 

The photon luminosity in model A is calculated based 
on the radiative flux mildly collimated within the po- 
lar angle of l9 = tan -^(25i ?s/60i?s) ^ 23° as we have 
mentioned in Section 14.2.11 Hence, the flows in model 
A would be identifled as extremely luminous objects of 
~ 13I/ph = 22Le for the accretion rate of ~ 300Le/c^, 
or 1.5 X 10^^ ergs~^, for the case of a black hole of 5OM0, 
in the face-on view, if an observer assumes the isotropic 
radiation fleld. Even higher apparent (isotropic) lumi- 
nosities are feasible for higher accretion rates, sin ce there 
is practically no limit to the mass accretion rate ( Ohsugal 
120061). 

The outflowing matter in model A is very hot, 
Tgas >10^ K, and dense, Tes ~ 1 with Tes being the Thom- 
son scattering optical depth, while the disk is relatively 
cold, Tgas - 10^-^ K (see Figures [5] and [g . Thus, the 
hot outflowing matter would Compton upscatt e r pho- 
tons from the disk surface. iKawashima et all ()2009l ) 
have shown that at high luminosities, Lph^I/E, a very 
hard spectral state appeared due to the effective inverse 
Compton scattering by the radiatively driven outflow. 
Thei r result can explain the hard X-ray spectra of ULXs 
(e.g..' Berghea et al.it200&) and seems to correspond to the 
ultraluminous state of ULXs (jGladstone et al.ir2009f ). In 
addition, since an inner part of the disk is obscured by 
the strong outflo w, the hot thermal comp onent might be 
invisible. Indeed. [Gladstone et al.l (|2009') have concluded 
based on the X-ray observations that the strong outflow 
envelopes the inner part of the supercritical disk, produc- 
ing the cool spectral component in some ULXs. Optically 
thick (res^3) low temperature (^ 10 keV) corona are ob- 
serv ed boTh in GRS 1915+10 5 and ULX, Holmberg IX 
X-1 (|Vierdavanti et al.l 120 lOal lbh. 

The geometrically thin, cold disk in model B is sur- 
rounded by the hot and rarefled matter (Figure[9]). Thus, 
the hot electrons in the corona would Compton upscat- 
ter the seed photons from the cold disk. While the 
thermal component of the spectra in active galactic nu- 
clei (AGNs) and BHBs is thought to be of disk origin, 
the Compton upscattering in the corona can account for 
the non-thermal hard component ( Deufel fc Sp ruit 20QQI: 
iLiu et al.l[2QQl iDone fc Gier lihski" '2004 l2005f ). No such 
extended disk corona was reproduced by the simulations 
of the local patch of the disk, which only show small, 
hot spo ts (with Tgas/^ 10 ^ K) appearing near the disk 
surface (jHirose et al.ll2006l ). 

In model B simulation, the optically thick, geometri- 
cally thin disk is truncated around r ~ 7i?s (see Section 
I4.3.ip . Such a disk truncat ion has been r eported by re- 
cent observations ( Kubota fc Pond [20041 : iTomsick et al.l 
I2OO9I : lYaniada et a l. 2009) and by the oretical spectral 
modeling (i Kawabata fc Mineshigelf20"TQ[ ). The blackbody 
radiation cannot be emitted within the truncation radius. 
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but instead the seed photons generated in the outer cool 
disk will be Compton upscattered by the hot electrons 
in the inner region. It then follows that the iron line is 
produced not at the innermost stable circular orbit but 
at the inner edge of the cool disk (truncation radius). If 
so, it would be diffi cult to derive the inforni ation about 
the black hole spin (jPone fc Gierlihskill2QQ6f ). 

Note, however, that evaporation of the disk gas due to 
the thermal conduction is not taken into account in the 
present simulation, although it could be an important 
ingredient to cause disk truncation (see Section 9.2.3 of 
[Kato et al .1120081 for comprehensive discussion). 

The disk in model C is optically thin and is composed 
of hot rarefied plasma. Synchrotron self-absorption, 
synchrotron self-Compton, and free — free eniission 
should be dominant r adiative processes (jLiu et al.ll2"QQl 
lOhsuga et aT '2005a"; "Kato et al.' "2009). Such spectra 
are observed in low-luminosity AGNs and in low— hard 
state of BHBs. Sgr A* is also a good example. Note 
that the electron temperature has been reported to de- 
viate from the ion t e mperature (iNaravan fc Yil 119951 : 
iNakamura et al.l [1996"; 'Manmoto et al."1997), although 
we assume one-temperature plasma in the present study. 
The emissio n of nonthermal electrons is thought to be 
non-trivial (|Yuan et al.l l2003f ). More detailed study is, 
however, beyond the scope of this paper. 

5.3. Future Work 

Global RMHD simulation is a relatively new research 
field, and, so there are plenty of issues to be explored 
in future. We stress again that higher resolution sim- 
ulations are needed to study the detailed structure of 
the geometrically thin disk like model B. The three- 
dimensional simulations should also be explored in future 
work. It is well known that no magnetic dynamo works 
in the axisymmetric calculations. Hence, three dimen- 
sional simulations are indispensable to make progress. 
We also stress that it is better to relax reflection sym- 
metry with respect to the equatorial plane, since such 
symmetry prevents the flows from going across the equa- 
torial plane. We plan such simulations as future work. In 
addition, we should take into account the relativistic ef- 
fects, since the outflow velocity becomes 10% — 70% of the 
speed of light near the rotation axis. In the present work, 
we start calculations with a rotating torus, in which the 
magnetic fields are purely poloidal. We should inves- 
tigate the inflow-outflow structure starting with other 
initial conditions in future. 

We used the FLD technique in the present study. This 
is quite a useful technique and gives appropriate radia- 
tion fields within the optically thick disk, but we should 
keep in mind that it does not always give precise ra- 
diation fields near the disk surface and above the disk 
(where the optical depth is around unity or less). In 
order to investigate the validity of the FLD, we calcu- 
late accurate radiative flux (RT flux) by solving radia- 
tion transfer equation along 10^ light rays at some points. 
The resulting RT flux is presented by a white vector in 
Figure [131 where the radiative flux evaluated based on 
the FLD method (FLD flux) is overlayed by a black vec- 
tor. Here we employ time- averaged (6—7 s) structure of 
model A, since the radiation force (radiative flux) in this 
model plays important roles. We find in this figure that 
the FLD flux shifts from the RT flux in the outflow re- 
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Fig. 13. — Radiative flux vectors by FLD approximation (black; 
FLD flux) and by integrating of radiation intensity (white; RT 
flux). Here, we solve the radiation transfer equation along 10^ 
light rays at each point for RT flux. Color contour indicates the 
radiation energy density, and the photosphere is plotted by dashed 
line. 

gions and near the photosphere, although the FLD flux 
is almost equal to RT flux in the optically thick disk re- 
gions. Especially, the radial component of the FLD flux 
is negative (inward flux), but, in contrast, that of the 
RT flux is positive (outward flux, e.g., see the fluxes at 
r ~ 7Rs with z ~ 12Rs). The radiation from the vicinity 
of the black hole works to produce the outward RT flux. 
However, the inward flux appears by the FLD method, 
since the direction of the flux is determined by the gra- 
dient of the radiation energy density. Thus, we should 
employ a more accurate method for the radiation fields in 
future. One of the i mproved methods is the so-called Ml 
closure scheme (Gonzalez e t al.l [20071 ). In this method, 
the 0th and 1st moment equations of radiation transfer 
are solved to update the radiation energy and the radia- 
tive flux. A closure relation, which is used to evaluate 
the radiation pressure tensor, is described in terms of 
the radiation energy and radiative flux (In the FLD ap- 
proximation, by contrast, the radiation pressure tensor 
is prescribed by the radiation energy only, and no infor- 
mation regarding radiative energy flux is used. This may 
result in somewhat inaccurate evaluation regarding the 
direction of radiation force in the jet acceleration region.) 

Throughout the present study, we use the frequency- 
integrated energy equation of the radiation. In or- 
der to obtain emergent spectra, we should perform the 
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frequency- dependent RMHD simulations. However, such 
simulations are so heavy that it is practically impossible 
to perform. Further, we only took into consideration the 
gas— radiation interaction via free— free and bound— free 
emission/absorption in our work. However, synchrotron 
emission/absorption and Compton scattering are impor- 
tant radiative processes. They would play an important 
role in the hot and tenuous regions found in accretion 
flows as disks with lower accretion rate (model C) and 
outflows in all models. 



6. CONCLUSIONS 

With a two-dimensional global RMHD code, we could 
reproduce three distinct inflow-outflow modes around 
black holes by adjusting a density normalization. Our 
three models with high, moderate, and low density nor- 
malizations correspond to the two-dimensional RMHD 
version of the slim disk (supercritical flow), the standard 
disk, and the RIAF, all with substantial outflows. 

We flnd the supercritical disk accretion flow, of which 
the photon luminosity exceeds the Eddington luminosity, 
in the case of the high density normalization (model A). 
The vertical component of the radiation force balances 
with that of the gravity in the disk region of z/r <2 but it 

largely exceeds the gravity above the disk, z/r>2. Our 

RMHD simulations reveals a new type of jet; i.e., the ra- 
diatively driven, magnetically collimated outflow, which 
might account for the jets of radio- loud NLSls and mi- 
croquasars. The disk, of which temperature is around 
10^~^ K, is surrounded by the hot outflowing matter, 
> 10^ K, which would induce the Compton upscatter- 
ing and the obscuration of the inner part of the disk. 
Because of the mildly collimated radiative flux, the ap- 
parent (isotropic) photon luminosity is 22Le, which is 
1.5 X 10^^ ergs~^ for the black hole of SOM©, in the face- 
on view. Even higher isotropic luminosity is feasible, if a 
greater amount of material is supplied to the black hole. 
Our supercritical model will be able to resolve the issue 
of the central engine of ULXs. 

When the moderate density normalization is employed, 
the radiative cooling is so effective that the cold geomet- 
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rically thin disk forms. This cold disk with ~ 10^ K 
is truncated at around 7i?s and enveloped by the hot 
rarefied atmosphere with Tgas > 10^ K, Compton up- 
scattering the seed photons from the cold disk. The cold 
thermal component and the non-thermal hard compo- 
nent of the spectra are observed in luminous AGNs and 
in the high-soft state of BHBs. The disk wind appears 
above and below the disk, which was not predicted in 
the framework of the standard disk model. The mag- 
netic pressure in the vertical direction is responsible for 
launching the gas from the disk surface. This result is 
consistent with the observations of the blueshifted ab- 
sorption lines. 

The simulations with low density normalization corre- 
sponds to the RIAF. The magnetic-pressure force in co- 
operation with the gas-pressure force drives the outflows. 
The flow releases the energy via jets rather than via ra- 
diation. The accretion flow as well as the outflows are 
hot and optically thin. Thus, the spectra would resemble 
those of the low-luminosity AGNs and of the BHBs their 
low-hard state. 

Finally, our simulations show that the vertically aver- 
aged disk visc osity is roug hly proportional to the total 
pressure (see lOhsuga et al.l l2QQ9). However, the local en- 
ergy dissipation rate is not simply proportional to the gas 
(and radiation) pressures in the vertical direction at least 
in the models with high and low density normalization. 
The vertical proflle of the dissipation rate is similar to 
that of the magnetic energy. This implies that the tradi- 
tional viscosity model is not perfectly correct. 
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